#-------------------------------------------------------------------------------
# Name:        Emin for IR/UV
# Purpose:
#
# Author:      gazizov
#
# Created:     29-11-2012
# Copyright:   (c) gazizov 2012
# Licence:     <your licence>
#-------------------------------------------------------------------------------

#=================================================================
from pylab import *

#=================================================================

def PlotBetaeepIR():
    beepIR = loadtxt('DataIR/beepIR.dat',unpack=True)
    Ep = exp(beepIR[0])
    b0 = beepIR[1]
    db0 = b0*beepIR[2]

    plot(Ep,b0,label=r'$\beta^{ee}(E)/H_0$')
    plot(Ep,db0,label=r'$H_0^{-1}\frac{d\beta^{ee}(E)}{d\ln(E)}$')
    xlabel('E, eV')
    xlim(1E14,1E23)
    ylabel(r'$\beta^{ee}(E)/H_0$, $H_0^{-1}d\beta^{ee}(E)/d\ln(E)$')
    xscale('log')
    text(1E21,0.011,'z=0',fontsize=16,color='r',backgroundcolor='w')
    txt = r'$p+\gamma_{\mathrm{EBL}} \rightarrow {e^+} + {e^-} + p$'
    text(4E19,0.01,txt,fontsize=16,color='g',backgroundcolor='w')
    legend(loc='upper left',ncol=1,shadow=False)
    grid(True)
    show()
    return

##PlotBetaeepIR()
##quit()


#=================================================================
##m_e       = 5.11E-4   # mass of electron in GeV
##m_p       = 0.939     # mass of nucleon in GeV
##ymax      = 1E1       # maximum energy of photon spectrum in eV
####Emineep   = 1E18*m_e*(m_e + m_p)/ymax  # in eV
##Egtheep = 1.0226E-3
##Emineep = 1E18*m_p*Egtheep/2E1
##lgEmineep = log10(Emineep)
##Egammin   = 0.15      # in GeV
##EminnX    = 1E18*m_p*Egammin/2E1    # in eV
##lgEminnX  = log10(EminnX)
##print lgEmineep
##print lgEminnX


def H(z):
    return sqrt(0.27*(1+z)**3 + 0.73)  # red-shift

def PlotPB():
    beepIR = loadtxt('DataIR/beepIR.dat',unpack=True)
    Ep = exp(beepIR[0])
    b0 = beepIR[1]
    Egg, Pp, Bp, Pn, Bn = loadtxt('Data/Check.dat',
    usecols=(0,1,2,4,5),unpack=True)
    plot(Ep,b0,linewidth=1.5)
    plot(Egg,Pp,label=r'$P_{pp}$')
    plot(Egg,Pn,label=r'$P_{pn}$')
    plot(Egg,Bp,label=r'$\beta_{pp}$')
    plot(Egg,Bn,label=r'$\beta_{pn}$')
    xscale('log')
    yscale('log')
    xlabel('E, eV')
    ylabel(r'$P(E)/H_0$, $\beta(E)/H_0$')
    xlim(1E15,1E23)
    ylim(1E-5,20.0)
    txteep = r'$p+\gamma_{\mathrm{EBL}} \rightarrow e^+e^-p$'
    text(2.4E19,1E-2,txteep,fontsize=15,
    rotation=-27,color='b') #,backgroundcolor='w')
    text(7.8E17,7.4,'z=0',fontsize=16,color='r'
    ,backgroundcolor='w')
    text(7E17,3.5E-5,'EBL by T. Kneiske et al.',
    fontsize=16,backgroundcolor='w')
    txtpg = r'$p+\gamma_{\mathrm{EBL}}\rightarrow p/n+X$'
    text(7.1E19,7.4,txtpg,fontsize=15,backgroundcolor='w')
    legend(loc='upper left',ncol=1,shadow=False)

    grid(True)
    show()
    return

def PlotBetaeep():
    for i in range(11):
        z   = 0.5*i
        fact = (1+z)**3/H(z)
        lbl = 'z='+str(z)
        if i == 0 or i == 10:
            lwd = 2
        else:
            lwd=1
        plot(Ep,fact*betaeeir[1+i*5],linewidth=lwd,label=lbl)
    xscale('log')
    yscale('log')
    xlabel('E, eV')
    ylabel(r'$\beta^{ee}_{\mathrm{EBL}}$ (E,z)')
    ylim(1E-11,10.0)
    txt1 = 'Relative energy losses on EBL in physical volume'
    text(3.E14,2.6,txt1,backgroundcolor='w',fontsize=13)
    txt2 = r'$p+\gamma_{\mathrm{EBL}}\rightarrow \, e^{+}+e^{-}+p$'
    text(4.2E13,9.2E-2,txt2,fontsize=15,color='b',backgroundcolor='w')
    txt3 = r'$\beta^{ee}_{\mathrm{EBL}}(E,z) = \frac{1}{H(z)}\, \frac{d\,\ln(E)}{dt}(E,z) $'
    text(9E15,5.7E-8,txt3,backgroundcolor='w',fontsize=16)
    legend(loc='lower center',ncol=3,shadow=False)

    grid(True)
    show()
    return


##PlotPB()
PlotBetaeep()